PROGRAM test6j

!  Test the RRF program, and specifically the 6j symbol routine, using the
!  sum rules in Brink & Satchler, Appendix II.

!  Suppress import of debug from rrfdata
USE rrf_module
USE wigner
USE input, ONLY: read_line, reada, readi, readu, reread, die, upcase
IMPLICIT NONE


INTEGER :: a, b, c, d, e, f, g, p, x, minj=0, maxj=12, test
TYPE(rrf), POINTER :: k=>null(), k1=>null(), k2=>null(), sum=>null(), &
    z=>null(), unity=>null()
LOGICAL :: eof, verbose=.false., printall, single=.false., allok, ok, &
    debug=.false.
CHARACTER(LEN=16) :: key
! CHARACTER(LEN=30) :: args
CHARACTER(LEN=8) :: verdict
CHARACTER(LEN=80) :: string(2)

write (0,"(a)") "Please wait ..."

call init_rrf(.true.,1000000)
print "(a,i0)", "Prime table includes primes up to ", prime(nprime)

! print "(a,i0)", "Long integers have kind ", long

do
  call read_line(eof)
  if (eof) exit
  call reada(key)
  select case(upcase(key))
    case default
      call die("Don't understand "//trim(key),.true.)
    case("EXIT","FINISH","QUIT")
      exit
    case("!","")
      cycle
    case("VERBOSE")
      verbose=.true.
    case("QUIET")
      verbose=.false.
    case ("PRINT")
      call readu(key)
      select case(key)
      case("ALL")
        printall=.true.
      case("FAILED","ERRORS")
        printall=.false.
      end select
    case("SINGLE")
      single=.true.
    case("MIN")
      call readi(minj)
      minj=2*minj
    case("MAX")
      call readi(maxj)
      maxj=2*maxj
    case("LIST","TABULATE")

      !  Tabulate
      !  a - f are double the actual j value
      do a=minj,maxj
        do b=a,maxj
          do d=0,maxj
            do e=abs(a-b),min(a+b,maxj),2
              do c=abs(d-e),min(d+e,maxj),2
                do f=abs(a-c),min(a+c,maxj),2
                  x=2
                  call sixj(a, b, c, d, e, f, x, k)
                  if ( nonzero(k) ) then
                    p=1
                    args=""
                    call put(a, args,p)
                    call put(b, args,p)
                    call put(c, args,p)
                    call put(d, args,p)
                    call put(e, args,p)
                    call put(f, args,p)
                    p=1
                    string(1)=""
                    call char4i(k, string(1), p, 80)
                    print "(a,2x,a)", args, trim(string(1))
                  end if
                end do
              end do
            end do
          end do
        end do
      end do

    case("SUM-RULE","SUMRULE","SUM", "TEST")
      !  SUM-RULE { 1 | 2 | 3 | 4 }
      !  If SINGLE has been specified, then
      !  SUM-RULE 4 [debug] a b c d f g
      !  where a-g are angular momentum values x 2
      test=1
      ! do while (item<nitems)
        call readu(key)
        select case(key)
        case("RULE")
          call readi(test)
        case default
          call reread(-1)
          call readi(test)
        end select
      ! end do
      select case(test)

      case(1)
        !  Sum rule 1
        !
        !  sum_c (-)^{2c} (2c+1) {a b c} = 1
        !                        {a b f}
        call unit(unity)
        x=2
        print "(/a)", "Sum rule 1"
        if (printall) print "(/a)", " a    b    f"
        allok=.true.
        !  a - f are double the actual j value
        do a=minj,maxj
          do b=a,maxj
            do f=abs(a-b),min(a+b,maxj),2
              !  Debug printing is normally off. If the LHS of the
              !  sum rule doesn't match the RHS, the calculation
              !  for that case is repeated with debug turned on.
              call sumrule1(a,b,f,x,ok,.false.)
              if (.not. ok) then
                call sumrule1(a,b,f,x,ok,.true.)
              end if
            end do
          end do
        end do
        if (allok) print "(a)", "Test completed -- no errors"

      case(2)
        !  Sum rule 2
        !
        !  sum_k (-)^{a+b+k} (2k+1) {a b k} = \delta_{f0){(2a+1)(2b+1)}^{1/2}
        !                           {b a f} 
        sum=>null()
        z=>null()
        x=2
        print "(/a)", "Sum rule 2"
        if (printall) print "(/a)", " a    b    f"
        !  a - f are double the actual j value
        do a=minj,maxj
          do b=a,maxj
            d=b
            e=a
            do f=abs(a-b),min(a+b,maxj),2
              call setzero(sum)
              do c=abs(a-b),abs(a+b),2
                call sixj(a, b, c, d, e, f, x, k)
                if (verbose) then
                  p=1
                  args=""
                  call put(a, args,p)
                  call put(b, args,p)
                  call put(c, args,p)
                  call put(d, args,p)
                  call put(e, args,p)
                  call put(f, args,p)
                  p=1
                  string(1)=""
                  call char4i(k, string(1), p, 80)
                  print "(a,2x,a)", args, trim(string(1))
                end if
                call int_to_rrf(int(c+1,long),z)
                call mult(k,z)
                if (mod((a+b+c)/2,2) .ne. 0) call chsign(k)
                call add(sum,k)
              end do
              if (f==0) then
                !  a & b are double the actual j value
                call int_to_rrf(int((a+1)*(b+1),long),z)
                call root(z)
              else
                call setzero(z)
              end if
              if (equal(sum,z)) then
                verdict="  O.K."
                ok=.true.
              else
                verdict=" Wrong!"
                ok=.false.
                if (allok .and. .not. printall) print "(/a)", " a    b    f"
                allok=.false.
              end if
              if (printall .or. .not. ok) then
                args=""
                p=1
                call put(a, args,p)
                call put(b, args,p)
                call put(f, args,p)
                p=1
                string(1)=""
                call char4i(sum, string(1), p, 80)
                print "(4a)", args(1:15), verdict, "Sum = ", trim(string(1))
              end if
            end do
          end do
        end do
        if (allok) print "(a)", "Test completed -- no errors"

      case(3)
        !  Sum rule 3
        !
        ! sum_e (2e+1)(2f+1) {a b e}{a b e} = \delta_{fg}
        !                    {c d f}{c d g}
        !
        print "(/a)", "Sum rule 3"
        if (printall) print "(/a)", " a    b    c    d    f    g"
        sum=>null()
        z=>null()
        call unit(unity)
        x=2
        allok=.true.
        !  a - g are double the actual j value
        do a=minj,maxj
          do b=a,maxj
            do c=0,maxj
              do d=0,maxj
                if (mod(c+d,2) .ne. mod(a+b,2)) cycle
                do f=max(abs(a-d),abs(b-c)),min(a+d,b+c,maxj),2
                  do g=f,min(a+d,b+c,maxj),2
                      !  Debug printing is normally off. If the LHS of the
                      !  sum rule doesn't match the RHS, the calculation
                      !  for that case is repeated with debug turned on.
                      call sumrule3(a,b,c,d,f,g,x,ok,.false.)
                      if (.not. ok) then
                        call sumrule3(a,b,c,d,f,g,x,ok,.true.)
                      end if
                  end do
                end do
              end do
            end do
          end do
        end do
        call clearall
        if (allok) print "(a)", "Test completed -- no errors"

      case(4)
        !  Sum rule 4
        !
        !  \sum_k (-)^{f+g+k} (2k+1) {a b k}{a b k} = {a d f}
        !                            {c d f}{d c g}   {b c g}
        !
        print "(/a)", "Sum rule 4"
        if (printall) print "(/a)", " a    b    c    d    f    g"
        call clear(sum)
        call clear(z)
        call unit(unity)
        x=2
        if (single) then
          call reada(key)
          print "(a)", key
          if (upcase(key) .eq. "DEBUG") then
            debug=.true.
          else
            debug=.false.
            call reread(-1)
          end if
          call readargs(a,b,c,d,f,g)
          call sumrule4(a,b,c,d,f,g,x,ok,debug)
          if (.not. ok .and..not. debug) call sumrule4(a,b,c,d,f,g,x,ok,.true.)
        else
          allok=.true.
          !  a - g are twice the actual j  value
          do a=minj,maxj
            do b=a,maxj
              do d=0,maxj
                do c=0,maxj
                  if (mod(c+d,2) .ne. mod(a+b,2)) cycle
                  do f=max(abs(a-d),abs(b-c)),min(a+d,b+c,maxj),2
                    do g=max(abs(a-c),abs(d-b)),min(a+c,d+b,maxj),2
                      !  Debug printing is normally off. If the LHS of the
                      !  sum rule doesn't match the RHS, the calculation
                      !  for that case is repeated with debug turned on.
                      call sumrule4(a,b,c,d,f,g,x,ok,.false.)
                      if (.not. ok) then
                        allok=.false.
                        call sumrule4(a,b,c,d,f,g,x,ok,.true.)
                      end if
                    end do
                  end do
                end do
              end do
            end do
          end do
        call clearall
        if (allok) print "(a)", "Test completed -- no errors"
      end if
    end select
  end select
  flush(6)
end do

CONTAINS

SUBROUTINE sumrule1(a, b, f, x, ok, debug)

!  Sum rule 1
!
!  sum_c (-)^{2c} (2c+1) {a b c} = 1
!                        {a b f}
!
!  Arguments a,b,f are each twice the angular momentum value.
!  x=2 always but is present for consistency with other routines.
INTEGER, INTENT(IN) :: a, b, f, x
LOGICAL, INTENT(OUT) :: ok
LOGICAL, INTENT(IN) :: debug
TYPE(rrf), POINTER :: k=>null(), sum=>null()
INTEGER :: c

! print "(a,6i3)", "sum rule 1: ", a,b,f
call setzero(sum)
do c=abs(a-b),abs(a+b),2
  call sixj(a, b, c, a, b, f, x, k)
  if (debug) then
    p=1
    args=""
    call put(a, args,p)
    call put(b, args,p)
    call put(c, args,p)
    call put(a, args,p)
    call put(b, args,p)
    call put(f, args,p)
    p=1
    string(1)=""
    call char4i(k, string(1), p, 80)
    print "(a,2x,a)", args, trim(string(1))
  end if
  call int_to_rrf(int(c+1,long),z)
  call mult(k,z)
  if (debug) then
    p=1
    string(1)=""
    call char4i(z, string(1), p, 80)
    p=p+3
    call char4i(k, string(1), p, 80)
    print "(a)", trim(string(1))
  end if
  call add(sum,k)
end do
!  (-)^{2k} = (-)^{a+b}
if (mod(a+b,2) .ne. 0) call chsign(sum)
if (equal(sum,unity)) then
  verdict="  O.K."
  ok=.true.
else
  verdict=" Wrong!"
  ok=.false.
  if (allok .and. .not. printall) print "(/a)", " a    b    f"
  allok=.false.
end if
if (printall .or. .not. ok) then
  p=1
  args=""
  call put(a, args, p)
  call put(b, args, p)
  call put(f, args, p)
  p=1
  string(1)=""
  call char4i(sum, string(1), p, 80)
  print "(4a)", args(1:15), verdict, "  Sum = ", trim(string(1))
end if

END SUBROUTINE sumrule1

!-----------------------------------------------------------------------

SUBROUTINE sumrule3(a, b, c, d, f, g, x, ok, debug)

!  Sum rule 3
!
! sum_e (2e+1)(2f+1) {a b e}{a b e} = \delta_{fg}
!                    {c d f}{c d g}
!
!  Arguments a,b,c,d,f,g are each twice the angular momentum value.
!  x=2 always but is present for consistency with other routines.
INTEGER, INTENT(IN) :: a, b, c, d, f, g, x
LOGICAL, INTENT(OUT) :: ok
LOGICAL, INTENT(IN) :: debug
TYPE(rrf), POINTER :: k=>null(), k1=>null(), k2=>null(), sum=>null()
INTEGER :: e
CHARACTER(LEN=6) :: buffer=""

args=""
p=1
call put(a, args,p)
call put(b, args,p)
call put(c, args,p)
call put(d, args,p)
call put(f, args,p)
call put(g, args,p)
! print "(a,6i3)", "sumrule 3: ", a,b,c,d,f,g
call setzero(sum)
do e=a+b,abs(a-b),-2
  call sixj(a, b, e, c, d, f, x, k1)
  call sixj(a, b, e, c, d, g, x, k2)
  call int_to_rrf(int(e+1,long),k)
  if (debug) then
    call put(e, buffer, p,.true.)
    print "(/2a)", "e = ", trim(buffer)
    call display(k,1,"k:")
    call display(k1,1,"k1:")
    call display(k2,1,"k2:")
  end if
  call mult(k,k1)
  call mult(k,k2)
  if (debug) then
    call display(k,1,"term:")
  end if
  call add(sum,k,debug)
  if (debug) then
    call display(sum,1,"sum:")
    if (sum%factor .ne. one) print "(a,i0)", "factor = ", sum%factor
  end if
end do
call multi(sum,int(f+1,long),1)
if (f .ne. g .and. iszero(sum) .or. f == g .and. equal(sum,unity)) then
  verdict="  O.K."
  ok=.true.
else
  verdict=" Wrong!"
  ok=.false.
  if (allok .and. .not. printall)                  &
      print "(/a)", " a    b    c    d    f    g"
  allok=.false.
end if
if (printall .or. .not. ok) then
  p=1
  string(1)=""
  call char4i(sum, string(1), p, 80)
  print "(4a)", args, verdict, "Sum = ", trim(string(1))
end if
if (debug) print "(1x)"

END SUBROUTINE sumrule3

!-----------------------------------------------------------------------

SUBROUTINE sumrule4(a, b, c, d, f, g, x, ok, debug)

!  Arguments a,b,c,d,f,g are each twice the angular momentum value.
!  x=2 always but is present for consistency with other routines.
INTEGER, INTENT(IN) :: a, b, c, d, f, g, x
LOGICAL, INTENT(OUT) :: ok
LOGICAL, INTENT(IN) :: debug
TYPE(rrf), POINTER :: k=>null(), k1=>null(), k2=>null(), sum=>null()
INTEGER :: e, p
CHARACTER(LEN=6) :: buffer=""

! print "(a,6i3)", "sumrule 4: ", a,b,c,d,f,g
call setzero(sum)
do e=max(abs(a-b),abs(c-d)),min(a+b,c+d),2   !  e <=> k in formula
  call sixj(a, b, e, c, d, f, x, k1)
  call sixj(a, b, e, d, c, g, x, k2)
  call int_to_rrf(int(e+1,long),k)
  if (debug) then
    call put(e, buffer, p,.true.)
    print "(/2a)", "e = ", trim(buffer)
    call display(k,1,"k:")
    call display(k1,1,"k1:")
    call display(k2,1,"k2:")
  end if
  call mult(k,k1)
  call mult(k,k2)
  if (mod((e+f+g)/2,2) .ne. 0) call chsign(k)
  if (debug) then
    call display(k,1,"term:")
  end if
  call add(sum,k,debug)
  if (debug) then
    call display(sum,1,"sum:")
    print "(a,i0)", "factor = ", sum%factor
  end if
end do
call sixj(a, d, f, b, c, g, x, z)
if (equal(sum,z)) then
  verdict="  O.K."
  ok=.true.
else
  verdict=" Wrong!"
  ok=.false.
  if (allok .and. .not. printall)                  &
      print "(/a)", " a    b    c    d    f    g"
  allok=.false.
end if
if (printall .or. .not. ok .or. single) then
  call put(a, args,p,.true.)
  call put(b, args,p)
  call put(c, args,p)
  call put(d, args,p)
  call put(f, args,p)
  call put(g, args,p)
  p=1
  string=""
  call char4i(sum, string(1), p, 80)
  p=1
  call char4i(z, string(2), p, 80)
  print "(a,t30,3a,t60,2a)", trim(args), verdict,               &
      "Sum = ", trim(string(1)), "6j = ", trim(string(2))
end if

END SUBROUTINE sumrule4

!-----------------------------------------------------------------------

SUBROUTINE readargs(a,b,c,d,e,f,g,h,i)

INTEGER, INTENT(OUT) :: a,b,c,d,e,f
INTEGER, INTENT(OUT), OPTIONAL :: g,h,i

call readarg(a)
call readarg(b)
call readarg(c)
call readarg(d)
call readarg(e)
call readarg(f)
if (present(g)) call readarg(g)
if (present(h)) call readarg(h)
if (present(i)) call readarg(i)

END SUBROUTINE readargs

!-----------------------------------------------------------------------

SUBROUTINE readarg(a)

INTEGER, INTENT(OUT) :: a
INTEGER :: p, x
CHARACTER(LEN=10) :: buff

call reada(buff)
p=1
x=2
call number(buff,p,a,x)
if (x==1) a=2*a

END SUBROUTINE readarg

SUBROUTINE put(j, s,p, reset)

CHARACTER(LEN=*) :: s
INTEGER, INTENT(IN) :: j
INTEGER, INTENT(INOUT) :: p
LOGICAL, INTENT(IN), OPTIONAL :: reset

if (present(reset)) then
  if (reset) then
    s=""
    p=1
  end if
end if

if (mod(j,2) == 0) then
  write (unit=s(p:p+1),fmt="(i2)") j/2
else
  write (unit=s(p:p+3),fmt="(i2,a2)") j,"/2"
endif
p=p+5

END SUBROUTINE put

!-----------------------------------------------------------------------

SUBROUTINE clearall

call clear(k)
call clear(k1)
call clear(k2)
call clear(sum)
call clear(z)
call clear(unity)

END SUBROUTINE clearall

END PROGRAM test6j
